1 Setup

1.1 Libs

library("vegan")
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-6.1
library("phyloseq")
library("ggplot2")
library("cowplot")
# #install.packages("colorBlindness")
# library("colorBlindness")
#BiocManager::install("microbiome")
#library("microbiome")
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
#install.packages("microViz", repos = c(davidbarnett = "https://david-barnett.r-universe.dev", getOption("repos")))
library(microViz)
## microViz version 0.12.4 - Copyright (C) 2021-2024 David Barnett
## ! Website: https://david-barnett.github.io/microViz
## ✔ Useful?  For citation details, run: `citation("microViz")`
## ✖ Silence? `suppressPackageStartupMessages(library(microViz))`
library("ggpubr")
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:microViz':
## 
##     stat_chull
## The following object is masked from 'package:cowplot':
## 
##     get_legend
# library('tidyverse')
library("microshades")
library("stringr")
library("speedyseq")
## 
## Attaching package: 'speedyseq'
## The following objects are masked from 'package:phyloseq':
## 
##     filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom,
##     tip_glom, transform_sample_counts
setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Mosquito_business/mosquito_microbes/04.microbiome_analysis/04b.comm_comp")

1.2 Data

ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.less.rds")
ps.clean #161 taxa, 458 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 161 taxa and 458 samples ]:
## sample_data() Sample Data:        [ 458 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 161 taxa by 11 taxonomic ranks ]:
## taxa are columns
ps.mq1 <- subset_samples(ps.clean,type=="A.albopictus")
ps.mq <- prune_taxa(taxa_sums(ps.mq1) > 0, ps.mq1)
ps.mq
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 89 taxa and 195 samples ]:
## sample_data() Sample Data:        [ 195 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 89 taxa by 11 taxonomic ranks ]:
## taxa are columns
ps.mw1 <- subset_samples(ps.clean,type=="Microbial Water")
ps.mw <- prune_taxa(taxa_sums(ps.mw1) > 0, ps.mw1)
ps.mw
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 49 taxa and 211 samples ]:
## sample_data() Sample Data:        [ 211 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 49 taxa by 11 taxonomic ranks ]:
## taxa are columns
##datasets but trimmed for mixed model things
ps.trim.mq <- readRDS("../04c.mixed_models/ps.trim.mq.rds")
ps.trim.mq #29 taxa, 195 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 29 taxa and 195 samples ]:
## sample_data() Sample Data:        [ 195 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 29 taxa by 11 taxonomic ranks ]:
## taxa are columns
ps.trim.mw <- readRDS("../04c.mixed_models/ps.trim.mw.rds")
ps.trim.mw #39 taxa, 211 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 39 taxa and 211 samples ]:
## sample_data() Sample Data:        [ 211 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 39 taxa by 11 taxonomic ranks ]:
## taxa are columns

2 PCOA

2.1 All samples

By ‘all’ here I really mean without separating the two experimental types (mosquitoes & mesocosm water)

2.1.1 Relative abundance & Bray

ps.clean.rel <- transform_sample_counts(ps.clean, function(x) x / sum(x))
ps.clean.rel
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 161 taxa and 458 samples ]:
## sample_data() Sample Data:        [ 458 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 161 taxa by 11 taxonomic ranks ]:
## taxa are columns
ps.all.rel1 <- subset_samples(ps.clean.rel,type=="A.albopictus"|type=="Microbial Water")
ps.all.rel1
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 161 taxa and 406 samples ]:
## sample_data() Sample Data:        [ 406 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 161 taxa by 11 taxonomic ranks ]:
## taxa are columns
ps.all.rel <- prune_taxa(taxa_sums(ps.all.rel1) > 0, ps.all.rel1)
ps.all.rel #103 taxa, 406 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 103 taxa and 406 samples ]:
## sample_data() Sample Data:        [ 406 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 103 taxa by 11 taxonomic ranks ]:
## taxa are columns
ord.all.rel <- ordinate(ps.all.rel, "PCoA", "bray")

plot.all.rel <- plot_ordination(ps.all.rel, ord.all.rel,color="type")+
  stat_ellipse()+
  theme_cowplot()
plot.all.rel #beautiful U shape

2.1.2 Aitchison (spelling?)

2.1.2.1 Full dataset

ps.clean.exp <- subset_samples(ps.clean,type=="A.albopictus"|type=="Microbial Water")

ps.exp.clr <- tax_transform(ps.clean.exp,"clr")

ord.exp.clr <- ordinate(ps.exp.clr, "PCoA", "euclidean")

gg.exp <- plot_ordination(ps.exp.clr, ord.exp.clr,color="type",label="mesocosm",shape="temperature")+
  scale_color_manual(name="Type",values=c("grey20","grey50"),labels=c("Mosq.","Water"))+
  theme_cowplot()+
  scale_shape_manual(name="Temperature",values=c(15,17),labels=c("Cool","Warm"))+
  ggtitle("Full dataset")
  #scale_fill_manual(values=c("white","black"),name="Type",labels=c("Mosq.","Water"))
gg.exp

##extract data so I can make it look how I want
df.exp.plot <- gg.exp[["data"]]
##new variable for plotting colors
df.exp.plot$inf.temp <- paste0(df.exp.plot$infusion,"_",df.exp.plot$temperature)
  
gg.exp.full <- ggplot(df.exp.plot,aes(x=Axis.1,y=Axis.2,fill=inf.temp,shape=inf.temp))+
  theme_cowplot()+
  stat_ellipse(aes(linetype=type,color=inf.temp),linewidth=0.6)+
  xlab("Axis 1 (52.3%)")+
  ylab("Axis 2 (11.4%)")+
  geom_point(color="gray70")+
  scale_color_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
  scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
  ggtitle("Both sample types")+
  scale_linetype_manual(values=c("solid","dashed"),name="Sample type",labels=c("Mosq.","Water"))+
  theme(legend.key.width=unit(1,"cm"),plot.title = element_text(face = "plain"))
  #theme(axis.text.x=element_text(angle=90),legend.key.width=unit(1,"cm"))
gg.exp.full

2.2 Mosquitoes - less trimmed

2.2.1 Multivariate plot

#ps.mq.rel <- transform_sample_counts(ps.trim.mq, function(x) x / sum(x) )
#ps.mq.rel <- transform_sample_counts(ps.mq, function(x) x / sum(x) )
#ps.trim.mq.hel <- transform_sample_counts(ps.trim.mq, function(OTU) sqrt(OTU/sum(OTU)))
ps.mq.clr <- tax_transform(ps.mq,"clr")

ord.mq.clr <- ordinate(ps.mq.clr, "PCoA", "euclidean")

gg.mq <- plot_ordination(ps.mq.clr, ord.mq.clr,color="infusion")+
  stat_ellipse()+
  theme_cowplot()
gg.mq

##extract data so I can make it look how I want
df.plot.mq <- gg.mq[["data"]]
##new variable for plotting colors
df.plot.mq$inf.temp <- paste0(df.plot.mq$infusion,"_",df.plot.mq$temperature)

gg.pcoa.aitch.mq <- ggplot(df.plot.mq,aes(x=Axis.1,y=Axis.2,fill=inf.temp,shape=inf.temp))+
  theme_cowplot()+
  stat_ellipse(aes(color=inf.temp),linewidth=0.6)+
  xlab("Axis 1 (17.6%)")+
  ylab("Axis 2 (14.3%)")+
  geom_point(color="gray70")+
  scale_color_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
  scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
  ggtitle("Mosquitoes only")+
  #scale_linetype_manual(values=c("solid","dashed","solid","dashed","solid","dashed"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  theme(legend.key.width=unit(1,"cm"),plot.title = element_text(face = "plain"))
  #theme(axis.text.x=element_text(angle=90),legend.key.width=unit(1,"cm"))

gg.pcoa.aitch.mq

2.2.2 Both panels - with diversity (Fig. 2)

Not knitting because requires the diversity script to be run also

#ggarrange(gg.exp.full,gg.pcoa.aitch.mq,labels=c("c","d"),common.legend=TRUE,legend="right")

gg.comp <- ggarrange(gg.exp.full,gg.pcoa.aitch.mq,labels=c("c","d"),common.legend=TRUE,legend="right")

ggarrange(gg.div,gg.comp,nrow=2,heights=c(0.8,1))

#ggsave("fig.div.alpha.beta.pdf",height=6,width=8)

2.3 Comparing prevalence-trimmed levels

ps.trim.exp <- merge_phyloseq(ps.trim.mq,ps.trim.mw)

ps.trim.exp.clr <- tax_transform(ps.trim.exp,"clr")

ord.trim.exp.clr <- ordinate(ps.trim.exp.clr, "PCoA", "euclidean")

gg.exp.trim <- plot_ordination(ps.trim.exp.clr, ord.trim.exp.clr,color="type",label="mesocosm",shape="temperature")+
  scale_color_manual(name="Type",values=c("grey20","grey50"),labels=c("Mosq.","Water"))+
  theme_cowplot()+
  scale_shape_manual(name="Temperature",values=c(15,17),labels=c("Cool","Warm"))+
  ggtitle("Prevalence-trimmed")
  #scale_fill_manual(values=c("white","black"),name="Type",labels=c("Mosq.","Water"))
gg.exp.trim

ggarrange(gg.exp,gg.exp.trim,common.legend=T,legend="right",labels=c("a","b"))

#ggsave(file="pcoa.mesoslabeled.png",width=9,height=4)
#ggsave(file="pcoa.mesoslabeled.pdf",width=9,height=4)

3 Trimmed vs. untrimmed stats

3.1 Mesocosm summed

3.1.1 Data - trimmed for mixed models

##merge by mesocosm - mosquitoes
ps.trim.mq.mer <- merge_samples(ps.trim.mq, "mesocosm")
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
samdf.trim.mq <- data.frame(sample_data(ps.trim.mq))
samdf.trim.mq.less <- samdf.trim.mq %>%
  dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.trim.mq.uniq <- distinct(samdf.trim.mq.less)
row.names(samdf.trim.mq.uniq) <- samdf.trim.mq.uniq$mesocosm
sample_data(ps.trim.mq.mer) <- sample_data(samdf.trim.mq.uniq)

##merge by mesocosm - water
ps.trim.mw.mer <- merge_samples(ps.trim.mw, "mesocosm")
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
samdf.trim.mw <- data.frame(sample_data(ps.trim.mw))
samdf.trim.mw.less <- samdf.trim.mw %>%
  dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.trim.mw.uniq <- distinct(samdf.trim.mw.less)
row.names(samdf.trim.mw.uniq) <- samdf.trim.mw.uniq$mesocosm
sample_data(ps.trim.mw.mer) <- sample_data(samdf.trim.mw.uniq)

3.1.2 Statz trimmed - mosquitoes

ps.trim.mq.mer.clr <- tax_transform(ps.trim.mq.mer,"clr")

seq.trim.mq.mer.clr <- data.frame(ps.trim.mq.mer.clr@otu_table)
samdf.trim.mq.mer.clr <- data.frame(ps.trim.mq.mer.clr@sam_data)

dist.aich.trim.mq.mer <- vegdist(seq.trim.mq.mer.clr, method="euclidean")

##beta - infusion
bet.aich.trim.mq.mer.inf <- betadisper(dist.aich.trim.mq.mer,samdf.trim.mq.mer.clr$infusion)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mq.mer.inf,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     2  10.345  5.1723 0.9967    999  0.372
## Residuals 55 285.407  5.1892
permutest(bet.aich.trim.mq.mer.inf,permutations=999,pairwise=T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     2  10.345  5.1723 0.9967    999  0.377
## Residuals 55 285.407  5.1892                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         OL      SG    SW
## OL         0.92400 0.151
## SG 0.92100         0.314
## SW 0.16761 0.33496
boxplot(bet.aich.trim.mq.mer.inf) #ns

##beta - temperature
bet.aich.trim.mq.mer.tem <- betadisper(dist.aich.trim.mq.mer,samdf.trim.mq.mer.clr$temperature)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mq.mer.tem,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     1   3.026  3.0261 0.5443    999  0.462
## Residuals 56 311.330  5.5595
boxplot(bet.aich.trim.mq.mer.tem) #ns

##beta - dispersal
bet.aich.trim.mq.mer.dis <- betadisper(dist.aich.trim.mq.mer,samdf.trim.mq.mer.clr$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mq.mer.dis,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     1   4.359  4.3591 0.8315    999   0.35
## Residuals 56 293.560  5.2421
boxplot(bet.aich.trim.mq.mer.dis) #ns

##adonis
adonis2(dist.aich.trim.mq.mer ~ infusion+temperature+dispersal, data=samdf.trim.mq.mer.clr, permutations=999) 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.aich.trim.mq.mer ~ infusion + temperature + dispersal, data = samdf.trim.mq.mer.clr, permutations = 999)
##             Df SumOfSqs      R2      F Pr(>F)    
## infusion     2    759.9 0.15531 5.1886  0.001 ***
## temperature  1    133.2 0.02723 1.8193  0.045 *  
## dispersal    1    118.7 0.02426 1.6213  0.063 .  
## Residual    53   3881.3 0.79320                  
## Total       57   4893.2 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#             Df SumOfSqs      R2      F Pr(>F)    
# infusion     2    759.9 0.15531 5.1886  0.001 ***
# temperature  1    133.2 0.02723 1.8193  0.044 *  
# dispersal    1    118.7 0.02426 1.6213  0.066 .  
# Residual    53   3881.3 0.79320                  
# Total       57   4893.2 1.00000                  

3.1.3 Statz trimmed - water

ps.trim.mw.mer.clr <- tax_transform(ps.trim.mw.mer,"clr")

seq.trim.mw.mer.clr <- data.frame(ps.trim.mw.mer.clr@otu_table)
samdf.trim.mw.mer.clr <- data.frame(ps.trim.mw.mer.clr@sam_data)

dist.aich.trim.mw.mer <- vegdist(seq.trim.mw.mer.clr, method="euclidean")

##beta - infusion
bet.aich.trim.mw.mer.inf <- betadisper(dist.aich.trim.mw.mer,samdf.trim.mw.mer.clr$infusion)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mw.mer.inf,permutations=999) #sig***
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
## Groups     2  73.567  36.784 22.286    999  0.001 ***
## Residuals 69 113.885   1.651                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permutest(bet.aich.trim.mw.mer.inf,permutations=999,pairwise=T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
## Groups     2  73.567  36.784 22.286    999  0.001 ***
## Residuals 69 113.885   1.651                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##            OL         SG    SW
## OL            1.0000e-03 0.013
## SG 4.1802e-04            0.001
## SW 1.1631e-02 3.8675e-08
boxplot(bet.aich.trim.mw.mer.inf) #all sig

##beta - temperature
bet.aich.trim.mw.mer.tem <- betadisper(dist.aich.trim.mw.mer,samdf.trim.mw.mer.clr$temperature)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mw.mer.tem,permutations=999) #0.057 .
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)  
## Groups     1  12.387 12.3865 3.8437    999  0.051 .
## Residuals 70 225.580  3.2226                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(bet.aich.trim.mw.mer.tem) #ns

##beta - dispersal
bet.aich.trim.mw.mer.dis <- betadisper(dist.aich.trim.mw.mer,samdf.trim.mw.mer.clr$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mw.mer.dis,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     1   8.585  8.5847 2.5233    999  0.127
## Residuals 70 238.149  3.4021
boxplot(bet.aich.trim.mw.mer.dis) #ns

##adonis
adonis2(dist.aich.trim.mw.mer ~ infusion+temperature+dispersal, data=samdf.trim.mw.mer.clr, permutations=999) 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.aich.trim.mw.mer ~ infusion + temperature + dispersal, data = samdf.trim.mw.mer.clr, permutations = 999)
##             Df SumOfSqs      R2       F Pr(>F)    
## infusion     2   5766.5 0.64480 67.6368  0.001 ***
## temperature  1    213.1 0.02382  4.9982  0.004 ** 
## dispersal    1    107.4 0.01201  2.5196  0.034 *  
## Residual    67   2856.1 0.31936                   
## Total       71   8943.1 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#             Df SumOfSqs      R2       F Pr(>F)    
# infusion     2   5766.5 0.64480 67.6368  0.001 ***
# temperature  1    213.1 0.02382  4.9982  0.002 ** 
# dispersal    1    107.4 0.01201  2.5196  0.053 .  
# Residual    67   2856.1 0.31936                   
# Total       71   8943.1 1.00000                   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

3.1.4 Data - not trimmed

##merge by mesocosm - mosquitoes
ps.mq.mer <- merge_samples(ps.mq, "mesocosm")
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
samdf.mq <- data.frame(sample_data(ps.mq))
samdf.mq.less <- samdf.mq %>%
  dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.mq.uniq <- distinct(samdf.mq.less)
row.names(samdf.mq.uniq) <- samdf.mq.uniq$mesocosm
sample_data(ps.mq.mer) <- sample_data(samdf.mq.uniq)

##merge by mesocosm - water
ps.mw.mer <- merge_samples(ps.mw, "mesocosm")
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
samdf.mw <- data.frame(sample_data(ps.mw))
samdf.mw.less <- samdf.mw %>%
  dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.mw.uniq <- distinct(samdf.mw.less)
row.names(samdf.mw.uniq) <- samdf.mw.uniq$mesocosm
sample_data(ps.mw.mer) <- sample_data(samdf.mw.uniq)

3.1.5 Statz not trimmed - mosquitoes

ps.mq.mer.clr <- tax_transform(ps.mq.mer,"clr")

seq.mq.mer.clr <- data.frame(ps.mq.mer.clr@otu_table)
samdf.mq.mer.clr <- data.frame(ps.mq.mer.clr@sam_data)

dist.aich.mq.mer <- vegdist(seq.mq.mer.clr, method="euclidean")

##beta - infusion
bet.aich.mq.mer.inf <- betadisper(dist.aich.mq.mer,samdf.mq.mer.clr$infusion)
#anova(bet.rad.uk2i)
permutest(bet.aich.mq.mer.inf,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     2  10.42  5.2081 0.6266    999  0.571
## Residuals 55 457.13  8.3114
permutest(bet.aich.mq.mer.inf,permutations=999,pairwise=T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     2  10.42  5.2081 0.6266    999  0.546
## Residuals 55 457.13  8.3114                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         OL      SG    SW
## OL         0.48000 0.574
## SG 0.48625         0.341
## SW 0.55470 0.32896
boxplot(bet.aich.mq.mer.inf) #ns

##beta - temperature
bet.aich.mq.mer.tem <- betadisper(dist.aich.mq.mer,samdf.mq.mer.clr$temperature)
#anova(bet.rad.uk2i)
permutest(bet.aich.mq.mer.tem,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq Mean Sq     F N.Perm Pr(>F)
## Groups     1   0.01  0.0080 0.001    999  0.978
## Residuals 56 438.96  7.8385
boxplot(bet.aich.mq.mer.tem) #ns

##beta - dispersal
bet.aich.mq.mer.dis <- betadisper(dist.aich.mq.mer,samdf.mq.mer.clr$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.aich.mq.mer.dis,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     1   0.44  0.4433 0.0583    999  0.801
## Residuals 56 425.79  7.6034
boxplot(bet.aich.mq.mer.dis) #ns

##adonis
adonis2(dist.aich.mq.mer ~ infusion+temperature+dispersal, data=samdf.mq.mer.clr, permutations=999) 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.aich.mq.mer ~ infusion + temperature + dispersal, data = samdf.mq.mer.clr, permutations = 999)
##             Df SumOfSqs      R2      F Pr(>F)    
## infusion     2    888.5 0.11387 3.6025  0.001 ***
## temperature  1    206.8 0.02651 1.6772  0.027 *  
## dispersal    1    171.6 0.02200 1.3918  0.113    
## Residual    53   6535.6 0.83763                  
## Total       57   7802.5 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#             Df SumOfSqs      R2      F Pr(>F)    
# infusion     2    759.9 0.15531 5.1886  0.001 ***
# temperature  1    133.2 0.02723 1.8193  0.044 *  
# dispersal    1    118.7 0.02426 1.6213  0.066 .  
# Residual    53   3881.3 0.79320                  
# Total       57   4893.2 1.00000                  

3.1.6 Statz not trimmed - water

ps.mw.mer.clr <- tax_transform(ps.mw.mer,"clr")

seq.mw.mer.clr <- data.frame(ps.mw.mer.clr@otu_table)
samdf.mw.mer.clr <- data.frame(ps.mw.mer.clr@sam_data)

dist.aich.mw.mer <- vegdist(seq.mw.mer.clr, method="euclidean")

##beta - infusion
bet.aich.mw.mer.inf <- betadisper(dist.aich.mw.mer,samdf.mw.mer.clr$infusion)
#anova(bet.rad.uk2i)
permutest(bet.aich.mw.mer.inf,permutations=999) #sig***
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq Mean Sq      F N.Perm Pr(>F)    
## Groups     2 81.949  40.974 32.543    999  0.001 ***
## Residuals 69 86.876   1.259                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permutest(bet.aich.mw.mer.inf,permutations=999,pairwise=T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq Mean Sq      F N.Perm Pr(>F)    
## Groups     2 81.949  40.974 32.543    999  0.001 ***
## Residuals 69 86.876   1.259                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##            OL         SG    SW
## OL            1.0000e-03 0.003
## SG 3.2943e-05            0.001
## SW 2.7862e-03 2.8516e-10
boxplot(bet.aich.mw.mer.inf) #all sig

##beta - temperature
bet.aich.mw.mer.tem <- betadisper(dist.aich.mw.mer,samdf.mw.mer.clr$temperature)
#anova(bet.rad.uk2i)
permutest(bet.aich.mw.mer.tem,permutations=999) #0.057 .
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)  
## Groups     1  13.948 13.9477 4.7498    999   0.03 *
## Residuals 70 205.555  2.9365                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(bet.aich.mw.mer.tem) #ns

##beta - dispersal
bet.aich.mw.mer.dis <- betadisper(dist.aich.mw.mer,samdf.mw.mer.clr$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.aich.mw.mer.dis,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)
## Groups     1   7.431  7.4314 2.3515    999  0.142
## Residuals 70 221.221  3.1603
boxplot(bet.aich.mw.mer.dis) #ns

##adonis
adonis2(dist.aich.mw.mer ~ infusion+temperature+dispersal, data=samdf.mw.mer.clr, permutations=999) 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.aich.mw.mer ~ infusion + temperature + dispersal, data = samdf.mw.mer.clr, permutations = 999)
##             Df SumOfSqs      R2       F Pr(>F)    
## infusion     2   6047.0 0.62704 62.3525  0.001 ***
## temperature  1    233.9 0.02426  4.8240  0.007 ** 
## dispersal    1    113.9 0.01181  2.3491  0.051 .  
## Residual    67   3248.9 0.33689                   
## Total       71   9643.7 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# adonis2(formula = dist.aich.mw.mer ~ infusion + temperature + dispersal, data = samdf.mw.mer.clr, permutations = 999)
#             Df SumOfSqs      R2       F Pr(>F)    
# infusion     2   6353.4 0.63122 63.3973  0.001 ***
# temperature  1    240.6 0.02391  4.8023  0.005 ** 
# dispersal    1    114.0 0.01133  2.2749  0.053 .  
# Residual    67   3357.2 0.33355                   
# Total       71  10065.2 1.00000                   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4 Unifrac things

4.1 Making the phylo tree info

Ran the following once & not doing it while knitting because it takes time

4.1.1 Subset needed taxa

Only need the taxa that end up in mosquitoes & mesocosm water samples

library("dada2")
library("DECIPHER")
library("phangorn")

setwd("~/Library/CloudStorage/GoogleDrive-nicolagk@hawaii.edu/My Drive/Mosquito_business/Mosquito_microbes/04.microbiome_analysis/04b.comm_comp")

ps.trim <- readRDS("../../02.process_asvs/ps.clean.trim.less.rds")
ps.trim #161 taxa, 458 samples

ps.exp <- subset_samples(ps.trim,type=="A.albopictus"|type=="Microbial Water")
ps.exp #406 samples
ps.exp.no0 <- prune_taxa(taxa_sums(ps.exp)!=0,ps.exp)
ps.exp.no0 #103 taxa

4.1.2 Alignment & tree things

Tutorial from dada2 author here

tax.exp <- data.frame(tax_table(ps.exp.no0))
head(tax.exp)

seqs <- getSequences(tax.exp$Sequence)
names(seqs) <- row.names(tax.exp) # This propagates to the tip labels of the tree
names(seqs)

alignment <- AlignSeqs(DNAStringSet(seqs), anchor=NA)
phang.align <- phyDat(as(alignment, "matrix"), type="DNA")
dm <- dist.ml(phang.align)
treeNJ <- NJ(dm) # Note, tip order != sequence order
fit = pml(treeNJ, data=phang.align)

## negative edges length changed to 0!
fitGTR <- update(fit, k=4, inv=0.2)
fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
                     rearrangement = "stochastic", control = pml.control(trace = 0))
detach("package:phangorn", unload=TRUE)
#saveRDS(fitGTR, file="phylo.fit.all.rds")

4.1.3 Add tree to phyloseq objects

ps.exp.no0@phy_tree <- phy_tree(fitGTR$tree)
ps.exp.no0

#saveRDS(ps.exp.no0,file="ps.exp.tree.rds")

4.2 Unifrac objects

Can skip above stuff

setwd("~/Library/CloudStorage/GoogleDrive-nicolagk@hawaii.edu/My Drive/Mosquito_business/mosquito_microbes/04.microbiome_analysis/04b.comm_comp")

ps.tree <- readRDS("ps.exp.tree.rds")
ps.tree
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 103 taxa and 406 samples ]:
## sample_data() Sample Data:        [ 406 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 103 taxa by 11 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 103 tips and 101 internal nodes ]:
## taxa are columns
fitGTR <- readRDS("phylo.fit.all.rds")

ps.trim.mq@phy_tree <- phy_tree(fitGTR$tree)
ps.trim.mq
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 29 taxa and 195 samples ]:
## sample_data() Sample Data:        [ 195 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 29 taxa by 11 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 103 tips and 101 internal nodes ]:
## taxa are columns
ps.trim.mw@phy_tree <- phy_tree(fitGTR$tree)
ps.trim.mw
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 39 taxa and 211 samples ]:
## sample_data() Sample Data:        [ 211 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 39 taxa by 11 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 103 tips and 101 internal nodes ]:
## taxa are columns

4.3 Mesocosm summed

4.3.1 Data - trimmed for mixed models

ps.trim.mq.rel <- transform_sample_counts(ps.trim.mq, function(x) x / sum(x))
ps.trim.mw.rel <- transform_sample_counts(ps.trim.mw, function(x) x / sum(x))

##merge by mesocosm - mosquitoes
ps.mq.mer <- merge_samples(ps.trim.mq.rel, "mesocosm")
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
samdf.mq <- data.frame(sample_data(ps.trim.mq.rel))
samdf.mq.less <- samdf.mq %>%
  dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.mq.uniq <- distinct(samdf.mq.less)
row.names(samdf.mq.uniq) <- samdf.mq.uniq$mesocosm
sample_data(ps.mq.mer) <- sample_data(samdf.mq.uniq)
ps.mq.mer
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 29 taxa and 58 samples ]:
## sample_data() Sample Data:        [ 58 samples by 4 sample variables ]:
## tax_table()   Taxonomy Table:     [ 29 taxa by 11 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 29 tips and 28 internal nodes ]:
## taxa are columns
##merge by mesocosm - water
ps.mw.mer <- merge_samples(ps.trim.mw.rel, "mesocosm")
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
samdf.mw <- data.frame(sample_data(ps.trim.mw.rel))
samdf.mw.less <- samdf.mw %>%
  dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.mw.uniq <- distinct(samdf.mw.less)
row.names(samdf.mw.uniq) <- samdf.mw.uniq$mesocosm
sample_data(ps.mw.mer) <- sample_data(samdf.mw.uniq)
ps.mw.mer
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 39 taxa and 72 samples ]:
## sample_data() Sample Data:        [ 72 samples by 4 sample variables ]:
## tax_table()   Taxonomy Table:     [ 39 taxa by 11 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 39 tips and 37 internal nodes ]:
## taxa are columns

4.4 Statz mesocosm summed - mosquitoes

# Transform counts to relative abundances
ps.mq.mer.rel <- transform_sample_counts(ps.mq.mer, function(x) x / sum(x))
ps.mq.mer.rel
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 29 taxa and 58 samples ]:
## sample_data() Sample Data:        [ 58 samples by 4 sample variables ]:
## tax_table()   Taxonomy Table:     [ 29 taxa by 11 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 29 tips and 28 internal nodes ]:
## taxa are columns
#seq.mq.mer.rel <- data.frame(ps.mq.mer.rel@otu_table)
sam.mq.mer.rel <- data.frame(ps.mq.mer.rel@sam_data)

dist.wuni.mq.rel <- phyloseq::distance(ps.mq.mer.rel, method = "wunifrac")

##beta - infusion
bet.wuni.mq.rel.inf <- betadisper(dist.wuni.mq.rel,sam.mq.mer.rel$infusion)
## Warning in betadisper(dist.wuni.mq.rel, sam.mq.mer.rel$infusion): some squared
## distances are negative and changed to zero
#anova(bet.rad.uk2i)
permutest(bet.wuni.mq.rel.inf,permutations=999) #sig **
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)   
## Groups     2 0.20598 0.102992 5.283    999  0.007 **
## Residuals 55 1.07223 0.019495                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permutest(bet.wuni.mq.rel.inf,permutations=999,pairwise=T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)   
## Groups     2 0.20598 0.102992 5.283    999  0.006 **
## Residuals 55 1.07223 0.019495                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##           OL        SG    SW
## OL           0.7700000 0.002
## SG 0.7319054           0.028
## SW 0.0097278 0.0479558
boxplot(bet.wuni.mq.rel.inf) #SW highest

##beta - temperature
bet.wuni.mq.rel.tem <- betadisper(dist.wuni.mq.rel,sam.mq.mer.rel$temperature)
#anova(bet.rad.uk2i)
permutest(bet.wuni.mq.rel.tem,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.07522 0.075224 2.7138    999  0.101
## Residuals 56 1.55230 0.027720
boxplot(bet.wuni.mq.rel.tem) #ns

##beta - dispersal
bet.wuni.mq.rel.dis <- betadisper(dist.wuni.mq.rel,sam.mq.mer.rel$dispersal)
## Warning in betadisper(dist.wuni.mq.rel, sam.mq.mer.rel$dispersal): some squared
## distances are negative and changed to zero
#anova(bet.rad.uk2i)
permutest(bet.wuni.mq.rel.dis,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.07251 0.072508 2.6287    999  0.106
## Residuals 56 1.54463 0.027583
boxplot(bet.wuni.mq.rel.dis) #ns

##adonis
adonis2(dist.wuni.mq.rel ~ infusion+temperature+dispersal, data=sam.mq.mer.rel, permutations=999) 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.wuni.mq.rel ~ infusion + temperature + dispersal, data = sam.mq.mer.rel, permutations = 999)
##             Df SumOfSqs      R2      F Pr(>F)    
## infusion     2  0.42030 0.21125 7.7307  0.001 ***
## temperature  1  0.08700 0.04373 3.2005  0.049 *  
## dispersal    1  0.04151 0.02086 1.5270  0.229    
## Residual    53  1.44075 0.72415                  
## Total       57  1.98956 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#             Df SumOfSqs      R2      F Pr(>F)    
# infusion     2  0.42030 0.21125 7.7307  0.001 ***
# temperature  1  0.08700 0.04373 3.2005  0.039 *  
# dispersal    1  0.04151 0.02086 1.5270  0.240    
# Residual    53  1.44075 0.72415                  
# Total       57  1.98956 1.00000                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4.5 Statz mesocosm summed - water

# Transform counts to relative abundances
ps.mw.mer.rel <- transform_sample_counts(ps.mw.mer, function(x) x / sum(x))
ps.mw.mer.rel
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 39 taxa and 72 samples ]:
## sample_data() Sample Data:        [ 72 samples by 4 sample variables ]:
## tax_table()   Taxonomy Table:     [ 39 taxa by 11 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 39 tips and 37 internal nodes ]:
## taxa are columns
#seq.mw.mer.rel <- data.frame(ps.mw.mer.rel@otu_table)
sam.mw.mer.rel <- data.frame(ps.mw.mer.rel@sam_data)

dist.wuni.mw.rel <- phyloseq::distance(ps.mw.mer.rel, method = "wunifrac")
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## Otu0049 -- in the phylogenetic tree in the data you provided.
##beta - infusion
bet.wuni.mw.rel.inf <- betadisper(dist.wuni.mw.rel,sam.mw.mer.rel$infusion)
#anova(bet.rad.uk2i)
permutest(bet.wuni.mw.rel.inf,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.000202 0.00010112 0.0716    999  0.933
## Residuals 69 0.097406 0.00141168
permutest(bet.wuni.mw.rel.inf,permutations=999,pairwise=T)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.000202 0.00010112 0.0716    999  0.944
## Residuals 69 0.097406 0.00141168                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         OL      SG    SW
## OL         0.85900 0.695
## SG 0.83196         0.889
## SW 0.70404 0.87936
boxplot(bet.wuni.mw.rel.inf) #ns

##beta - temperature
bet.wuni.mw.rel.tem <- betadisper(dist.wuni.mw.rel,sam.mw.mer.rel$temperature)
#anova(bet.rad.uk2i)
permutest(bet.wuni.mw.rel.tem,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.00373 0.0037270 0.7637    999  0.379
## Residuals 70 0.34161 0.0048801
boxplot(bet.wuni.mw.rel.tem) #ns

##beta - dispersal
bet.wuni.mw.rel.dis <- betadisper(dist.wuni.mw.rel,sam.mw.mer.rel$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.wuni.mw.rel.dis,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.00071 0.0007117 0.1412    999  0.694
## Residuals 70 0.35294 0.0050421
boxplot(bet.wuni.mw.rel.dis) #ns

##adonis
adonis2(dist.wuni.mw.rel ~ infusion+temperature+dispersal, data=sam.mw.mer.rel, permutations=999) 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.wuni.mw.rel ~ infusion + temperature + dispersal, data = sam.mw.mer.rel, permutations = 999)
##             Df SumOfSqs      R2        F Pr(>F)    
## infusion     2  1.04644 0.77043 168.1578  0.001 ***
## temperature  1  0.10069 0.07413  32.3606  0.001 ***
## dispersal    1  0.00265 0.00195   0.8519  0.409    
## Residual    67  0.20847 0.15348                    
## Total       71  1.35825 1.00000                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#             Df SumOfSqs      R2        F Pr(>F)    
# infusion     2   3.7520 0.78212 168.4687  0.001 ***
# temperature  1   0.2880 0.06004  25.8669  0.001 ***
# dispersal    1   0.0111 0.00231   0.9955  0.330    
# Residual    67   0.7461 0.15552                    
# Total       71   4.7972 1.00000                    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5 Rel. abun. bar plot - regular~

Need to make some edits to ‘culture_name’ column for plotting. Specifically, adding wolbachia strain names

ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.less.rds")
ps.clean #161 taxa, 458 samples
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 161 taxa and 458 samples ]:
## sample_data() Sample Data:        [ 458 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 161 taxa by 11 taxonomic ranks ]:
## taxa are columns
tax.clean <- data.frame(ps.clean@tax_table)

#write.csv(tax.clean,file="tax.clean.toedit.csv")

##read in with edits
tax.edits <- read.csv("tax.clean.edits.csv",row.names=1)

ps.clean.copy <- ps.clean
ps.clean.copy@tax_table <- tax_table(as.matrix(tax.edits))

ps.mq1 <- subset_samples(ps.clean.copy,type=="A.albopictus")
ps.mq <- prune_taxa(taxa_sums(ps.mq1) > 0, ps.mq1)
ps.mq
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 89 taxa and 195 samples ]:
## sample_data() Sample Data:        [ 195 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 89 taxa by 11 taxonomic ranks ]:
## taxa are columns
ps.mw1 <- subset_samples(ps.clean.copy,type=="Microbial Water")
ps.mw <- prune_taxa(taxa_sums(ps.mw1) > 0, ps.mw1)
ps.mw
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 49 taxa and 211 samples ]:
## sample_data() Sample Data:        [ 211 samples by 24 sample variables ]:
## tax_table()   Taxonomy Table:     [ 49 taxa by 11 taxonomic ranks ]:
## taxa are columns

5.1 Prep - mosquitoes

##Having OTU & size columns in front of Kingdom messes up phyloseq stuff
tax.mq <- data.frame(ps.mq@tax_table)
tax.mq.cut <- tax.mq[,3:10]
tax.mq.cut$OTU <- row.names(tax.mq.cut)
##custom grouping for microshades below
tax.mq.cut$gen_culture <- paste0(tax.mq.cut$Genus,"_",tax.mq.cut$Pool_name)
##making a copy of ps object so I don't overwrite original
ps.mq.order <- ps.mq
ps.mq.order@tax_table <- tax_table(as.matrix(tax.mq.cut))

##relative abundance
ps.mq.rel <- transform_sample_counts(ps.mq.order, function(x) x / sum(x))

samdf.mq <- data.frame(ps.mq.rel@sam_data)
samdf.mq$glom <- paste0(samdf.mq$sex,"_",samdf.mq$infusion,"_",samdf.mq$temperature)
samdf.mq$glom
##   [1] "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C" "M_OL_C" "M_OL_C" "M_OL_C" "F_OL_C"
##   [9] "M_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C"
##  [17] "M_OL_C" "M_OL_H" "M_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "M_OL_H" "M_OL_H"
##  [25] "F_OL_H" "F_OL_H" "M_OL_H" "F_OL_H" "M_OL_H" "M_OL_C" "F_OL_C" "M_OL_C"
##  [33] "M_OL_C" "F_OL_C" "M_OL_C" "F_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C"
##  [41] "M_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_H" "M_OL_H" "F_OL_H" "F_OL_H"
##  [49] "F_OL_H" "F_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "F_OL_H" "F_OL_H" "M_OL_H"
##  [57] "M_OL_H" "F_OL_H" "M_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "F_OL_H" "M_OL_H"
##  [65] "M_OL_H" "F_OL_H" "F_OL_H" "M_OL_C" "F_OL_C" "M_OL_C" "F_OL_C" "F_OL_C"
##  [73] "F_OL_C" "M_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C" "M_OL_C" "M_OL_C"
##  [81] "M_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "F_OL_C" "F_OL_C" "F_OL_H" "M_OL_H"
##  [89] "M_OL_H" "F_OL_H" "F_OL_H" "M_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "M_OL_H"
##  [97] "F_OL_H" "F_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "M_SG_C" "M_SG_C" "M_SG_H"
## [105] "F_SG_C" "M_SG_C" "F_SG_C" "F_SG_C" "F_SG_C" "F_SG_C" "F_SG_C" "F_SG_C"
## [113] "M_SG_C" "M_SG_C" "M_SG_C" "M_SG_C" "F_SG_C" "F_SG_C" "M_SG_C" "F_SG_C"
## [121] "F_SG_H" "M_SG_H" "M_SG_H" "M_SG_H" "F_SG_H" "F_SG_H" "M_SG_H" "F_SG_H"
## [129] "F_SG_H" "F_SG_H" "M_SG_H" "M_SG_H" "F_SG_H" "F_SG_H" "F_SG_H" "M_SG_H"
## [137] "F_SG_H" "F_SG_H" "M_SG_H" "M_SG_C" "F_SG_H" "F_SG_H" "F_SG_H" "M_SG_H"
## [145] "M_SG_H" "F_SG_H" "F_SW_C" "M_SW_C" "F_SW_C" "M_SW_C" "M_SW_C" "M_SW_C"
## [153] "F_SW_C" "M_SW_C" "M_SW_C" "F_SW_H" "F_SW_H" "F_SW_H" "F_SW_C" "F_SW_C"
## [161] "M_SW_C" "M_SW_C" "F_SW_C" "M_SW_C" "M_SW_C" "F_SW_C" "F_SW_C" "F_SW_C"
## [169] "M_SW_C" "M_SW_C" "M_SW_C" "M_SW_C" "F_SW_C" "M_SW_C" "M_SW_C" "M_SW_C"
## [177] "F_SW_H" "M_SW_H" "F_SW_H" "M_SW_C" "M_SW_C" "M_SW_C" "M_SW_C" "F_SW_C"
## [185] "M_SW_C" "M_SW_C" "F_SW_C" "F_SW_C" "M_SW_C" "F_SW_H" "F_SW_H" "F_SW_H"
## [193] "F_SW_H" "M_SW_H" "F_SW_H"
ps.mq.rel@sam_data <- sample_data(samdf.mq)

ps.mq.rel.glom <- merge_samples2(ps.mq.rel, "glom")
ps.mq.rel.glom.rel <- transform_sample_counts(ps.mq.rel.glom, function(x) x / sum(x))

plot_bar(ps.mq.rel.glom.rel,fill="Family")
## Warning in psmelt(physeq, as = "data.frame"): The rank names: 
## OTU
##  have been renamed to: 
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.

##Who are the top families for microshades
# Merges ASVs that have the same taxonomy rank (Genus)
ps.mq.rel.fam <- tax_glom(ps.mq.rel, taxrank = "Family") #28 families, 39 gen.

# Calculate taxa sum
top5 = head(sort(colSums(otu_table(ps.mq.rel.fam)), decreasing = TRUE), 10)
# Combine count and taxonomyTable
top5 = cbind(as.data.frame(tax_table(ps.mq.rel.fam)[names(top5),]), Count = top5)
top5
##                       Class            Order             Family Genus Species
## Otu0001 Alphaproteobacteria    Rickettsiales    Anaplasmataceae  <NA>    <NA>
## Otu0003 Gammaproteobacteria  Xanthomonadales   Xanthomonadaceae  <NA>    <NA>
## Otu0005         Bacteroidia Flavobacteriales      Weeksellaceae  <NA>    <NA>
## Otu0002 Gammaproteobacteria Enterobacterales Enterobacteriaceae  <NA>    <NA>
## Otu0022             Bacilli  Lactobacillales  Carnobacteriaceae  <NA>    <NA>
## Otu0010 Gammaproteobacteria  Pseudomonadales   Pseudomonadaceae  <NA>    <NA>
## Otu0031 Alphaproteobacteria      Rhizobiales  Xanthobacteraceae  <NA>    <NA>
## Otu0121 Gammaproteobacteria Enterobacterales           Orbaceae  <NA>    <NA>
## Otu0007 Alphaproteobacteria  Acetobacterales   Acetobacteraceae  <NA>    <NA>
## Otu0137 Gammaproteobacteria  Burkholderiales   Burkholderiaceae  <NA>    <NA>
##         Sequence  OTU Pool_name gen_culture        Count
## Otu0001     <NA> <NA>      <NA>        <NA> 170.19659314
## Otu0003     <NA> <NA>      <NA>        <NA>  11.05234660
## Otu0005     <NA> <NA>      <NA>        <NA>   5.79876664
## Otu0002     <NA> <NA>      <NA>        <NA>   3.93453195
## Otu0022     <NA> <NA>      <NA>        <NA>   1.85376844
## Otu0010     <NA> <NA>      <NA>        <NA>   1.83582760
## Otu0031     <NA> <NA>      <NA>        <NA>   0.14403788
## Otu0121     <NA> <NA>      <NA>        <NA>   0.03718247
## Otu0007     <NA> <NA>      <NA>        <NA>   0.03185857
## Otu0137     <NA> <NA>      <NA>        <NA>   0.02061469
#Anaplasmataceae
#Xanthomonadaceae
#Weeksellaceae
#Enterobacteriaceae

5.2 Microshades - mosquitoes

mdf.mq <- prep_mdf(ps.mq.rel.glom.rel,subgroup_level="gen_culture")
## Warning in speedyseq::psmelt(.): The rank names: 
## OTU
##  have been renamed to: 
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.
# Create a color object for the specified data
col.mdf.mq <- create_color_dfs(mdf.mq, top_orientation = FALSE,group_level="Family",subgroup_level="gen_culture",selected_groups=c("Anaplasmataceae","Xanthomonadaceae","Weeksellaceae","Enterobacteriaceae"),cvd=TRUE)

#Extract
col.mdf.mq.m <- col.mdf.mq$mdf
col.mdf.mq.c <- col.mdf.mq$cdf

##default plot
plot_microshades(col.mdf.mq.m, col.mdf.mq.c)+
  #facet_wrap(scales="free")+
  theme_cowplot()+
  facet_wrap(~infusion,scales="free")+
  scale_x_discrete(labels=c("F_cool","F_warm","M_cool","M_warm"))+
  theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))+
  xlab("")+
  ggtitle("Mosquitoes")+
  guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))

# 
# ##sample data got removed during conglomerating
# col.mdf.mq.m$sex <- substr(col.mdf.mq.m$Sample,1,1)
# col.mdf.mq.m$infusion <- substr(col.mdf.mq.m$Sample,3,4)
# col.mdf.mq.m$temperature <- substr(col.mdf.mq.m$Sample,6,6)
# col.mdf.mq.m$infusion <- sub("SW","PW",col.mdf.mq.m$infusion)
# col.mdf.mq.m$infusion <- factor(col.mdf.mq.m$infusion,levels=c("OL","SG","PW"))

col.mdf.mq.c$group
##  [1] "Enterobacteriaceae-Kosakonia_KOSA1"                  
##  [2] "Enterobacteriaceae-Family_Enterobacteriaceae_ENTE1/2"
##  [3] "Enterobacteriaceae-Klebsiella_KLEB1"                 
##  [4] "Enterobacteriaceae-Family_Enterobacteriaceae_"       
##  [5] "Enterobacteriaceae-Other"                            
##  [6] "Weeksellaceae-Chryseobacterium_CHRY1"                
##  [7] "Weeksellaceae-Chryseobacterium_"                     
##  [8] "Xanthomonadaceae-Stenotrophomonas_STEN1"             
##  [9] "Xanthomonadaceae-Stenotrophomonas_STEN2"             
## [10] "Xanthomonadaceae-Stenotrophomonas_"                  
## [11] "Xanthomonadaceae-Thermomonas_XANT1"                  
## [12] "Anaplasmataceae-Wolbachia_wAlbB"                     
## [13] "Anaplasmataceae-Wolbachia_wAlbA"                     
## [14] "Anaplasmataceae-Wolbachia_"                          
## [15] "Other-Carnobacterium_CARN1"                          
## [16] "Other-Pseudomonas_PSEU2"                             
## [17] "Other-Bradyrhizobium_"                               
## [18] "Other-Asaia_ASAI1/2"                                 
## [19] "Other-Other"
##can't see these ones:
bye <- c("Enterobacteriaceae-Family_Enterobacteriaceae_", "Enterobacteriaceae-Other", "Weeksellaceae-Chryseobacterium_", "Xanthomonadaceae-Stenotrophomonas_", "Xanthomonadaceae-Thermomonas_XANT1", "Anaplasmataceae-Wolbachia_", "Other-Bradyrhizobium_", "Other-Asaia_ASAI1/2", "Other-Other")
col.mdf.mq.m.less <- col.mdf.mq.m[!col.mdf.mq.m$group %in% bye,]
##fix the lengths of these ones:
# col.mdf.mq.m.less$group <- sub("Enterobacteriaceae-Enterobacteriaceae_unclassified_ENTE1/2","Enterobacteriaceae-unclassified_ENTE1/2",col.mdf.mq.m.less$group)

##desired order
col.mdf.mq.m.less$group <- factor(col.mdf.mq.m.less$group,levels=c("Other-Carnobacterium_CARN1","Other-Pseudomonas_PSEU2","Anaplasmataceae-Wolbachia_wAlbB","Anaplasmataceae-Wolbachia_wAlbA","Xanthomonadaceae-Stenotrophomonas_STEN1","Xanthomonadaceae-Stenotrophomonas_STEN2","Weeksellaceae-Chryseobacterium_CHRY1","Enterobacteriaceae-Kosakonia_KOSA1","Enterobacteriaceae-Family_Enterobacteriaceae_ENTE1/2","Enterobacteriaceae-Klebsiella_KLEB1"))

##color order
col.mdf.mq.c$hex
##  [1] "#7D3560" "#A1527F" "#CC79A7" "#E794C1" "#EFB6D6" "#148F77" "#009E73"
##  [8] "#098BD9" "#56B4E9" "#7DCCFF" "#BCE1FF" "#9D654C" "#C17754" "#F09163"
## [15] "#616161" "#8B8B8B" "#B7B7B7" "#D6D6D6" "#F5F5F5"
hex.colors <- c("#616161","#8B8B8B","#9D654C","#C17754","#098BD9","#56B4E9","#148F77","#7D3560","#A1527F","#CC79A7")

gg.bars.mosq <- ggplot(data=col.mdf.mq.m.less,aes(x=Sample,y=Abundance,fill=group,color=group))+
  geom_bar(stat="identity", position="stack")+
  facet_wrap(~infusion,scales="free")+
  scale_fill_manual(values=hex.colors)+
  scale_color_manual(values=hex.colors)+
  ggtitle("Mosquitoes")+
  guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))+
  theme_cowplot()+
  scale_x_discrete(labels=c("F_cool","F_warm","M_cool","M_warm"))+
  xlab("")+
  theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))
gg.bars.mosq

5.3 Prep - water

##Having OTU & size columns in front of Kingdom messes up phyloseq stuff
tax.mw <- data.frame(ps.mw@tax_table)
tax.mw.cut <- tax.mw[,3:10]
tax.mw.cut$OTU <- row.names(tax.mw.cut)
##custom grouping for microshades below
tax.mw.cut$gen_culture <- paste0(tax.mw.cut$Genus,"_",tax.mw.cut$Pool_name)
##making a copy of ps object so I don't overwrite original
ps.mw.order <- ps.mw
ps.mw.order@tax_table <- tax_table(as.matrix(tax.mw.cut))

##relative abundance
ps.mw.rel <- transform_sample_counts(ps.mw.order, function(x) x / sum(x))

##new sample type so I can group them
samdf.mw <- data.frame(ps.mw.rel@sam_data)
samdf.mw$glom <- paste0(samdf.mw$day,"_",samdf.mw$infusion,"_",samdf.mw$temperature)
samdf.mw$glom
##   [1] "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_H" "12_OL_H" "12_OL_H"
##   [8] "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_H" "12_OL_H" "12_OL_H"
##  [15] "12_OL_H" "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_H" "12_OL_H"
##  [22] "12_OL_H" "12_OL_H" "12_SG_C" "12_SG_C" "12_SG_C" "12_SG_C" "12_SG_H"
##  [29] "12_SG_H" "12_SG_H" "12_SG_H" "12_SG_C" "12_SG_C" "12_SG_C" "12_SG_C"
##  [36] "12_SG_H" "12_SG_H" "12_SG_H" "12_SG_C" "12_SG_C" "12_SG_C" "12_SG_C"
##  [43] "12_SG_H" "12_SG_H" "12_SG_H" "12_SG_H" "12_SW_C" "12_SW_C" "12_SW_C"
##  [50] "12_SW_C" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_C" "12_SW_C"
##  [57] "12_SW_C" "12_SW_C" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_C"
##  [64] "12_SW_C" "12_SW_C" "12_SW_C" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_H"
##  [71] "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_H" "20_OL_H" "20_OL_H"
##  [78] "20_OL_H" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_H" "20_OL_H"
##  [85] "20_OL_H" "20_OL_H" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_H"
##  [92] "20_OL_H" "20_OL_H" "20_OL_H" "20_SG_C" "20_SG_C" "20_SG_C" "20_SG_C"
##  [99] "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_C" "20_SG_C" "20_SG_C"
## [106] "20_SG_C" "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_C" "20_SG_C"
## [113] "20_SG_C" "20_SG_C" "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_H" "20_SW_C"
## [120] "20_SW_C" "20_SW_C" "20_SW_C" "20_SW_H" "20_SW_H" "20_SW_H" "20_SW_C"
## [127] "20_SW_C" "20_SW_C" "20_SW_C" "20_SW_H" "20_SW_H" "20_SW_H" "20_SW_H"
## [134] "20_SW_C" "20_SW_C" "20_SW_C" "20_SW_C" "20_SW_H" "20_SW_H" "20_SW_H"
## [141] "20_SW_H" "4_OL_C"  "4_OL_C"  "4_OL_C"  "4_OL_C"  "4_OL_H"  "4_OL_H" 
## [148] "4_OL_H"  "4_OL_H"  "4_OL_C"  "4_OL_C"  "4_OL_C"  "4_OL_C"  "4_OL_H" 
## [155] "4_OL_H"  "4_OL_H"  "4_OL_H"  "4_OL_C"  "4_OL_C"  "4_OL_C"  "4_OL_C" 
## [162] "4_OL_H"  "4_OL_H"  "4_OL_H"  "4_OL_H"  "4_SG_C"  "4_SG_C"  "4_SG_C" 
## [169] "4_SG_C"  "4_SG_H"  "4_SG_H"  "4_SG_H"  "4_SG_C"  "4_SG_C"  "4_SG_C" 
## [176] "4_SG_H"  "4_SG_H"  "4_SG_H"  "4_SG_H"  "4_SG_C"  "4_SG_C"  "4_SG_C" 
## [183] "4_SG_C"  "4_SG_H"  "4_SG_H"  "4_SG_H"  "4_SG_H"  "4_SW_C"  "4_SW_C" 
## [190] "4_SW_C"  "4_SW_C"  "4_SW_H"  "4_SW_H"  "4_SW_H"  "4_SW_H"  "4_SW_C" 
## [197] "4_SW_C"  "4_SW_C"  "4_SW_C"  "4_SW_H"  "4_SW_H"  "4_SW_H"  "4_SW_H" 
## [204] "4_SW_C"  "4_SW_C"  "4_SW_C"  "4_SW_C"  "4_SW_H"  "4_SW_H"  "4_SW_H" 
## [211] "4_SW_H"
ps.mw.rel@sam_data <- sample_data(samdf.mw)

ps.mw.rel.glom <- merge_samples2(ps.mw.rel, "glom")
ps.mw.rel.glom.rel <- transform_sample_counts(ps.mw.rel.glom, function(x) x / sum(x))

plot_bar(ps.mw.rel.glom.rel,fill="Family")
## Warning in psmelt(physeq, as = "data.frame"): The rank names: 
## OTU
##  have been renamed to: 
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.

##Who are the top families for microshades
# Merges ASVs that have the same taxonomy rank 
ps.mw.rel.fam <- tax_glom(ps.mw.rel, taxrank = "Family") #28 families, 39 gen.

# Calculate taxa sum
top5 = head(sort(colSums(otu_table(ps.mw.rel.fam)), decreasing = TRUE), 10)
# Combine count and taxonomyTable
top5 = cbind(as.data.frame(tax_table(ps.mw.rel.fam)[names(top5),]), Count = top5)
top5
##                       Class            Order                 Family Genus
## Otu0002 Gammaproteobacteria Enterobacterales     Enterobacteriaceae  <NA>
## Otu0003 Gammaproteobacteria  Xanthomonadales       Xanthomonadaceae  <NA>
## Otu0007 Alphaproteobacteria  Acetobacterales       Acetobacteraceae  <NA>
## Otu0005         Bacteroidia Flavobacteriales          Weeksellaceae  <NA>
## Otu0010 Gammaproteobacteria  Pseudomonadales       Pseudomonadaceae  <NA>
## Otu0014 Gammaproteobacteria Enterobacterales            Erwiniaceae  <NA>
## Otu0031 Alphaproteobacteria      Rhizobiales      Xanthobacteraceae  <NA>
## Otu0001 Alphaproteobacteria    Rickettsiales        Anaplasmataceae  <NA>
## Otu0070 Gammaproteobacteria Enterobacterales Order_Enterobacterales  <NA>
## Otu0162 Gammaproteobacteria Enterobacterales         Morganellaceae  <NA>
##         Species Sequence  OTU Pool_name gen_culture       Count
## Otu0002    <NA>     <NA> <NA>      <NA>        <NA> 149.4822380
## Otu0003    <NA>     <NA> <NA>      <NA>        <NA>  32.0710953
## Otu0007    <NA>     <NA> <NA>      <NA>        <NA>  13.1670955
## Otu0005    <NA>     <NA> <NA>      <NA>        <NA>   6.4043253
## Otu0010    <NA>     <NA> <NA>      <NA>        <NA>   4.3610484
## Otu0014    <NA>     <NA> <NA>      <NA>        <NA>   3.3758270
## Otu0031    <NA>     <NA> <NA>      <NA>        <NA>   1.0753763
## Otu0001    <NA>     <NA> <NA>      <NA>        <NA>   0.7593486
## Otu0070    <NA>     <NA> <NA>      <NA>        <NA>   0.2676740
## Otu0162    <NA>     <NA> <NA>      <NA>        <NA>   0.0255795
#Enterobacteriaceae
#Xanthomonadaceae
#Acetobacteraceae
#Weeksellaceae

5.4 Microshades - water

mdf.mw <- prep_mdf(ps.mw.rel.glom.rel,subgroup_level="gen_culture")
## Warning in speedyseq::psmelt(.): The rank names: 
## OTU
##  have been renamed to: 
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.
# Create a color object for the specified data
col.mdf.mw <- create_color_dfs(mdf.mw, top_orientation = FALSE,group_level="Family",subgroup_level="gen_culture",selected_groups=c("Enterobacteriaceae","Xanthomonadaceae","Acetobacteraceae","Weeksellaceae"),cvd=TRUE)

#Extract
col.mdf.mw.m <- col.mdf.mw$mdf
col.mdf.mw.c <- col.mdf.mw$cdf

##default plot
plot_microshades(col.mdf.mw.m, col.mdf.mw.c)+
  #facet_wrap(scales="free")+
  theme_cowplot()+
  #scale_x_discrete(labels=c("F_cool","F_warm","M_cool","M_warm"))+
  theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))+
  xlab("")+
  guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))

col.mdf.mw.c.new <- color_reassign(col.mdf.mw.c,
                          group_assignment = c("Enterobacteriaceae","Xanthomonadaceae","Acetobacteraceae","Weeksellaceae"),
                          color_assignment = c("micro_cvd_purple","micro_cvd_blue","micro_cvd_green","micro_cvd_turquoise"),group_level="Family")

##check new colors
plot_microshades(col.mdf.mw.m, col.mdf.mw.c.new)+
  #facet_wrap(scales="free")+
  theme_cowplot()+
  #scale_x_discrete(labels=c("F_cool","F_warm","M_cool","M_warm"))+
  theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))+
  xlab("")+
  guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))

# ##sample data got removed during conglomerating
# col.mdf.mw.m2 <- cbind(col.mdf.mw.m,data.frame(str_split_fixed(col.mdf.mw.m$Sample,"_",3)))
# col.mdf.mw.m2$X1==col.mdf.mw.m2$day #day lines up
# #X2 = infusion, X3=temp

col.mdf.mw.m$infusion <- sub("SW","PW",col.mdf.mw.m$infusion)
col.mdf.mw.m$infusion <- factor(col.mdf.mw.m$infusion,levels=c("OL","SG","PW"))

col.mdf.mw.c.new$group
##  [1] "Weeksellaceae-Chryseobacterium_CHRY1"                
##  [2] "Acetobacteraceae-Asaia_ASAI1/2"                      
##  [3] "Acetobacteraceae-Asaia_ASAI3"                        
##  [4] "Acetobacteraceae-Asaia_"                             
##  [5] "Xanthomonadaceae-Stenotrophomonas_STEN1"             
##  [6] "Xanthomonadaceae-Stenotrophomonas_STEN2"             
##  [7] "Xanthomonadaceae-Stenotrophomonas_"                  
##  [8] "Xanthomonadaceae-Thermomonas_XANT1"                  
##  [9] "Enterobacteriaceae-Kosakonia_KOSA1"                  
## [10] "Enterobacteriaceae-Family_Enterobacteriaceae_ENTE1/2"
## [11] "Enterobacteriaceae-Kosakonia_KOSA2"                  
## [12] "Enterobacteriaceae-Klebsiella_KLEB1"                 
## [13] "Enterobacteriaceae-Other"                            
## [14] "Other-Pantoea_PANT1/2/3"                             
## [15] "Other-Pseudomonas_PSEU2"                             
## [16] "Other-Bradyrhizobium_"                               
## [17] "Other-Pseudomonas_"                                  
## [18] "Other-Other"
##can't see these ones:
bye.mw <- c("Acetobacteraceae-Asaia_","Xanthomonadaceae-Stenotrophomonas_","Xanthomonadaceae-Thermomonas_XANT1","Other-Pseudomonas_","Other-Other")
  
col.mdf.mw.m.less <- col.mdf.mw.m[!col.mdf.mw.m$group %in% bye.mw,]

# ##fix the lengths of these ones:
# ##Entero group
# col.mdf.mw.m.less$group <- sub("Enterobacteriaceae-Enterobacteriaceae_unclassified_ENTE1/2","Enterobacteriaceae-unclassified_ENTE1/2",col.mdf.mw.m.less$group)
# col.mdf.mw.m.less$group <- sub("Enterobacteriaceae-Enterobacteriaceae_unclassified_KLEB1","Enterobacteriaceae-Klebsiella_KLEB1",col.mdf.mw.m.less$group)
# col.mdf.mw.m.less$group <- sub("Enterobacteriaceae-Enterobacteriaceae_unclassified_KOSA2","Enterobacteriaceae-Kosakonia_KOSA2",col.mdf.mw.m.less$group)
##Other group
col.mdf.mw.m.less$group <- sub("Other-Bradyrhizobium_","Other-Bradyrhizobium",col.mdf.mw.m.less$group)
# col.mdf.mw.m.less$group <- sub("Other-Gammaproteobacteria_unclassified_","Other-Gammaproteobacteria_unclassified",col.mdf.mw.m.less$group)

##desired order
col.mdf.mw.m.less$group <- factor(col.mdf.mw.m.less$group,levels=c("Other-Pantoea_PANT1/2/3","Other-Pseudomonas_PSEU2","Other-Bradyrhizobium","Acetobacteraceae-Asaia_ASAI1/2","Acetobacteraceae-Asaia_ASAI3","Xanthomonadaceae-Stenotrophomonas_STEN1","Xanthomonadaceae-Stenotrophomonas_STEN2","Weeksellaceae-Chryseobacterium_CHRY1","Enterobacteriaceae-Kosakonia_KOSA1","Enterobacteriaceae-Family_Enterobacteriaceae_ENTE1/2","Enterobacteriaceae-Klebsiella_KLEB1","Enterobacteriaceae-Kosakonia_KOSA2","Enterobacteriaceae-Other"))

##switching some colors to match mosquitoes
hex.colors.mw <- c("#616161","#8B8B8B","#B7B7B7","#4E7705","#6D9F06","#098BD9","#56B4E9","#148F77","#7D3560","#A1527F","#CC79A7","#E794C1","#EFB6D6")

##Sample order
col.mdf.mw.m.less$Sample <- factor(col.mdf.mw.m.less$Sample, levels=c("4_OL_C","4_OL_H","12_OL_C","12_OL_H","20_OL_C","20_OL_H","4_SG_C","4_SG_H","12_SG_C","12_SG_H","20_SG_C","20_SG_H","4_SW_C","4_SW_H","12_SW_C","12_SW_H","20_SW_C","20_SW_H"))

gg.bars.water <- ggplot(data=col.mdf.mw.m.less,aes(x=Sample,y=Abundance,fill=group,color=group))+
  geom_bar(stat="identity", position="stack")+
  facet_wrap(~infusion,scales="free")+
  scale_fill_manual(values=hex.colors.mw)+
  scale_color_manual(values=hex.colors.mw)+
  ggtitle("Mesocosm water")+
  guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))+
  theme_cowplot()+
  scale_x_discrete(labels=c("D4_cool","D4_warm","D12_cool","D12_warm","D20_cool","D20_warm"))+
  xlab("")+
  theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))
gg.bars.water

5.5 Both mosq & water

ggarrange(gg.bars.mosq,gg.bars.water,labels=c("(a)","(b)"),nrow=2)

#ggsave("total_comp.pdf",width=10.5,height=7)
#ggsave("total_comp.png",width=10.5,height=7)

6 Rel. abun. box plots - mosquitoes

ps.trim.mq.rel <- transform_sample_counts(ps.trim.mq, function(x) x / sum(x))

melt.trim.mq.rel <- psmelt(ps.trim.mq.rel)
## Warning in psmelt(ps.trim.mq.rel): The rank names: 
## OTU
##  have been renamed to: 
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.

6.1 Wolb in sex

melt.trim.mq.rel.wolb <- subset(melt.trim.mq.rel,Genus=="Wolbachia")

ggplot(data = melt.trim.mq.rel.wolb, aes(x = sex, y = Abundance))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(aes(color = OTU), height = 0, width = .2) +
  facet_wrap(~OTU,scales="free")

melt.trim.mq.rel.1.4 <- subset(melt.trim.mq.rel, OTU %in% c("Otu0001", "Otu0004"))

melt.trim.mq.rel.1.4$Longer_name <- melt.trim.mq.rel.1.4$OTU
melt.trim.mq.rel.1.4$Longer_name <- sub("Otu0001","Otu0001_wAlbB",melt.trim.mq.rel.1.4$Longer_name)
melt.trim.mq.rel.1.4$Longer_name <- sub("Otu0004","Otu0004_wAlbA",melt.trim.mq.rel.1.4$Longer_name)

gg.sex.wolb <- ggplot(data = melt.trim.mq.rel.1.4, aes(x = sex, y = Abundance))+
  geom_jitter(width=0.2,alpha=0.5)+
  geom_boxplot(outlier.shape = NA,alpha=0.5)+
  facet_wrap(~Longer_name,scales="free")+
  xlab("Mosquito sex")+
  ylab("Relative abundance/sample")+
  theme_bw()+
  scale_x_discrete(labels=c("Female","Male"))+
  guides(color="none")
gg.sex.wolb

# ggsave("relabun.wolb.pdf",width=4,height=3)
# ggsave("relabun.wolb.png",width=4,height=3)

6.2 Chry & pseu in infusion

melt.trim.mq.rel.5.10 <- subset(melt.trim.mq.rel, OTU %in% c("Otu0005", "Otu0010"))

gg.inf <- ggplot(data = melt.trim.mq.rel.5.10, aes(x = infusion, y = Abundance))+
  geom_jitter(width=0.2,alpha=0.5)+
  geom_boxplot(outlier.shape = NA,alpha=0.5)+
  facet_wrap(~Longer_name,scales="free")+
  ylab("Relative abundance/sample")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  theme_bw()+
  guides(color="none")
gg.inf

#ggsave("relabun.chrypseu.pdf",width=4,height=2.5)

6.3 Kleb in temp

melt.trim.mq.rel.kleb <- subset(melt.trim.mq.rel, OTU %in% c("Otu0009","Otu0022"))

gg.tem <- ggplot(data = melt.trim.mq.rel.kleb, aes(x = temperature, y = Abundance))+
  geom_jitter(width=0.2,alpha=0.5)+
  geom_boxplot(outlier.shape = NA,alpha=0.5)+
  facet_wrap(~Longer_name,scales="free")+
  ylab("Relative abundance/sample")+
  theme_bw()+
  xlab("Temperature")+
  scale_x_discrete(labels=c("Cool","Warm"))+
  ylim(0,0.017)+
  guides(color="none")
gg.tem 
## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 182 rows containing missing values or values outside the scale range
## (`geom_point()`).

#1 in cool for kleb, 1 in warm for carn, 5 in cool for carn

#ggsave("relabun.chrypseu.pdf",width=4,height=2.5)

6.4 Multiple arranged

#ggarrange(gg.sex.wolb,gg.inf,gg.tem,nrow=3,labels=c("(a)","(b)","(c)"))
#ggsave(file="relabun.bars.all.pdf",width=5,height=8)

7 Relative abundance box plots for mesocosm water

ps.trim.mw.rel <- transform_sample_counts(ps.trim.mw, function(x) x / sum(x))

melt.trim.mw.rel <- psmelt(ps.trim.mw.rel)
## Warning in psmelt(ps.trim.mw.rel): The rank names: 
## OTU
##  have been renamed to: 
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.

7.1 Infusion differences

melt.trim.mw.inf <- subset(melt.trim.mw.rel, OTU %in% c("Otu0003","Otu0005","Otu0007","Otu0024"))

ggplot(data = melt.trim.mw.inf, aes(x = infusion, y = Abundance))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(aes(color = Longer_name), height = 0, width = .2) +
  facet_wrap(~Longer_name,scales="free")

gg.mw.inf <- ggplot(data = melt.trim.mw.inf, aes(x = infusion, y = Abundance))+
  geom_jitter(width=0.2,alpha=0.5)+
  geom_boxplot(outlier.shape = NA,alpha=0.5)+
  facet_wrap(~Longer_name,scales="free")+
  xlab("Infusion")+
  ylab("Relative abundance/sample")+
  theme_bw()+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  guides(color="none")
gg.mw.inf

#ggsave("relabun.wolb.pdf",width=4,height=3)

8 Session info

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.2.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Pacific/Honolulu
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] speedyseq_0.5.3.9021 stringr_1.5.1        microshades_1.13    
##  [4] ggpubr_0.6.0         microViz_0.12.4      dplyr_1.1.4         
##  [7] cowplot_1.1.3        ggplot2_3.5.1        phyloseq_1.48.0     
## [10] vegan_2.6-6.1        lattice_0.22-6       permute_0.9-7       
## 
## loaded via a namespace (and not attached):
##  [1] ade4_1.7-22             tidyselect_1.2.1        farver_2.1.2           
##  [4] Biostrings_2.72.1       fastmap_1.2.0           digest_0.6.37          
##  [7] lifecycle_1.0.4         cluster_2.1.6           survival_3.7-0         
## [10] magrittr_2.0.3          compiler_4.4.1          rlang_1.1.4            
## [13] sass_0.4.9              tools_4.4.1             igraph_2.0.3           
## [16] utf8_1.2.4              yaml_2.3.10             data.table_1.15.4      
## [19] knitr_1.48              ggsignif_0.6.4          labeling_0.4.3         
## [22] plyr_1.8.9              abind_1.4-5             Rtsne_0.17             
## [25] withr_3.0.1             purrr_1.0.2             BiocGenerics_0.50.0    
## [28] grid_4.4.1              stats4_4.4.1            fansi_1.0.6            
## [31] multtest_2.60.0         biomformat_1.32.0       colorspace_2.1-1       
## [34] Rhdf5lib_1.26.0         scales_1.3.0            iterators_1.0.14       
## [37] MASS_7.3-61             cli_3.6.3               rmarkdown_2.28         
## [40] crayon_1.5.3            generics_0.1.3          microbiome_1.26.0      
## [43] httr_1.4.7              reshape2_1.4.4          ape_5.8                
## [46] cachem_1.1.0            rhdf5_2.48.0            zlibbioc_1.50.0        
## [49] splines_4.4.1           parallel_4.4.1          XVector_0.44.0         
## [52] vctrs_0.6.5             Matrix_1.7-0            carData_3.0-5          
## [55] jsonlite_1.8.8          car_3.1-2               IRanges_2.38.1         
## [58] S4Vectors_0.42.1        rstatix_0.7.2           foreach_1.5.2          
## [61] tidyr_1.3.1             jquerylib_0.1.4         glue_1.7.0             
## [64] codetools_0.2-20        stringi_1.8.4           gtable_0.3.5           
## [67] GenomeInfoDb_1.40.1     UCSC.utils_1.0.0        munsell_0.5.1          
## [70] tibble_3.2.1            pillar_1.9.0            htmltools_0.5.8.1      
## [73] rhdf5filters_1.16.0     GenomeInfoDbData_1.2.12 R6_2.5.1               
## [76] evaluate_0.24.0         Biobase_2.64.0          highr_0.11             
## [79] backports_1.5.0         broom_1.0.6             bslib_0.8.0            
## [82] Rcpp_1.0.13             gridExtra_2.3           nlme_3.1-166           
## [85] mgcv_1.9-1              xfun_0.47               forcats_1.0.0          
## [88] pkgconfig_2.0.3